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Oh 

q Previous high resolution cosmological simulations predict the first 

£j stars to appear in the early universe to be very massive and to form 

in isolation. Here we discuss a cosmological simulation in which the 
central 5OM clump breaks up into two cores, having a mass ratio of 
i— h two to one, with one fragment collapsing to densities of 10 8 g cm 3 . 

The second fragment, at a distance of ~ 800 astronomical units, 
is also optically thick to its own cooling radiation from molecular 
hydrogen lines, but is still able to cool via collision-induced emission. 
The two dense peaks will continue to accrete from the surrounding 
cold gas reservoir over a period of ~ 10 5 years and will likely form a 
0> binary star system. 

^s! Hydrodynamical simulations that start from cosmological initial conditions have pre- 

dieted that the first luminous objects in the universe were isolated stars with masses in 
the range 30 — 3OOM , based on accretion rates calculated in the absence of accretion- 
inhibiting factors (1-3). Idealized protostellar evolution simulations suggest accretion 
should end when the star reaches ~ 100M o (4), but even though most relevant feedback 
mechanisms of the protostars on their accretion flow are known (5), no fully self-consistent 
radiation hydrodynamical simulations reaching the main sequence have yet been possible, 
even in one dimension. The early stages of proto-stellar evolution, however, are under- 
stood in some detail from spherically symmetric calculations (6, 7) as well as semi-analytic 
models (8). 

The environment in which the first stars form is calculated by following early-universe 
evolution of the primordial gas and dark matter; by this means, the first stars are found to 
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form in halos with a total mass of ~ 1O 6 M , with collapse driven first by cooling via molec- 
ular hydrogen ro- vibrational lines and later by collision induced emission (6,9, 10). While 
several tens of calculations (1,11, 12) have followed the collapse of a primordial protostel- 
lar halo to relatively low densities (~ 10 _12 g cm -3 ), only two previous calculations have 
followed the collapse from cosmological initial conditions to protostellar densities (2, 13). 
Both of these calculations utilized a chemical model that included the relevant medium- 
and high-density physics, such as heating from the formation of molecular hydrogen, col- 
lisionally induced emission, three-body molecular hydrogen formation, and gas opacity 
at high densities. Several simulations have followed the collapse of a single metal-free 
cloud starting from idealized initial conditions (3, 14, 15). Parameter studies of rotating 
cylinders have shown that metal- free gas can fragment (16) and have found fragmenta- 
tion of spherical distributions of extremely low- and zero-metallicity gas (17), including 
in 3D nested grid parameter studies (18). Additionally, studies of zero-metallicity gas 
collapsing in isolation have found fragmentation at low densities (~ 10" 16 g cm" 3 ) (19). 
However, no previous simulation starting from cosmological initial conditions has shown 
fragmentation in a cosmological context. 

We present simulations performed with the adaptive mesh refinement code Enzo 
(O'Shea et al.(2004))- These simulations were initialized at a redshift of 99 in a box 
300 kpc h" 1 (comoving). (Detailed simulation parameters and methods can be found 
in the supporting online material.) We stopped the simulation at a maximum (proper) 
baryon density of 1.61 x 10" 8 g cm" 3 at z — 19.08 (189 million years after the Big Bang) 
and 1.72 R Q (proper) peak spatial resolution, where we have resolved the first massive 
halo that forms in this simulation volume. The simulation presented here is explicitly a 
simulation of the first generation of Population III star formation. The second generation 
of Population III stars are still primordial in composition but are affected by previous 
generation of stars, potentially by effects such as kinetic energy injection, radiation back- 
grounds or cosmic rays. We have simulated five realizations of first-generation Population 
III stars that evolve to at least this density and have found fragmentation in one. The 
dark matter halo merger history and the evolution of baryonic properties at redshifts 
above 20 are typical of previous simulations (1,20). 

The collapsing halo fragmented at a density of ~ 2.0 x 10" 13 g cm" 3 at redshift 
z = 19.08, roughly. We have conducted resolution studies, where we found fragmen- 
tation but a change in the separation of the identified cores. (More information can 
be found in the supporting online materials.) To study the fragmentation of the cloud, 
we identified gravitationally-bound, topologically-connected regions at or above a single 
density within the halo's virial radius. The density at which these regions diverged was 
~ 2.0 x 10" 13 g cm -3 , however, because the difference was only a thin strip of lower density 
gas, we conducted all subsequent analysis on connected isodensity surfaces with minimum 
density of 3.0 x 10" 13 g cm -3 . We selected material within two isodensity surfaces, here- 
after referred to as Core A and Core B (Fig. 1). Core A is the more massive core, with 
a mass of 10.0M o , and Core B of 6.3M . Determining whether these two objects will 
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ultimately merge is nontrivial, but estimates can be made based on the current orbital 
mechanics as well as collapse conditions of each core in isolation. If these two objects form 
hydrostatically-supported protostellar cores and most of the intervening gas is accreted 
onto them, they will remain separate and form a binary star system. 

Both cores are oblate and are embedded in a crescent-shaped cloud. Even though they 
can not be reliably identified more than ~ 175 years before the end of the simulation, the 
beginning of fragmentation is already evident in both density and molecular hydrogen, 
where lower-density, mostly-atomic regions are permeated by filaments of higher-molecular 
fraction gas. (Fig. 1, top row.) Averaging over the time period during which the two cores 
are identifiable and distinct, Core A experiences 0.061 M /year mass flux across its density 
isosurface, while Core B has mass flux of 0.049 M /year. (Fig. 2) The flow across both 
density isosurfaces is roughly linear. The peak densities in both cores evolve at close to 
free-fall trajectories and Core B collapses ~ 60 years after Core A. 

The separation of the centers of mass of the two cores is 800AU. Within a minimally 
enclosing sphere of radius 780AU located at the center of mass of the two cores, we find 
a total mass of 52M ; within a sphere of radius twice the separation of the cores, we find 
a mass of 99 M . By defining 



we calculated the characteristic radius of each core. Core B, the lower-mass core, is at this 
time less dense then Core A, and has a smaller characteristic radius. Given its current 
angular momentum and assuming both angular momentum conservation and no material 
between the cores its final Keplerian radius is given by 



where Lb is the specific angular momentum of Core B in the rest frame of Core A 
and Ma is the mass of Core A. R^ cp evaluates to 2400 AU, which is three times their 
current distance, indicating that the velocity of Core B is super-Keplerian. The tangential 
velocity of Core B with respect to Core A is 5.9 km s _1 , and it has a radial velocity of 
—5.3 km s _1 . Given the current velocities, separation and masses of the cores and the 
assumption that the two cores exist with their current trajectories in a vacuum, a simple 
orbit integration shows that their closest approach of 480 AU would occur after 360 years, 
more than 10 dynamical times of the lighter core. However, this vacuum assumption could 
be misleading. Another limit can be derived assuming the remaining 46M mass between 
Core A and B instantaneously were accreted onto only Core A, then the orbit integrations 
show it would take over 400 years for B to come within 300 AU of A. Because the time 
scale of the formation of hydrostatically-supported, protostellar cores in either core is close 
to the free fall time, an accreting binary system will form. Following the formulation of 
the dynamical friction force on a massive perturber in a near-circular orbit (21), the 
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loss of angular momentum through dynamical friction does not substantially increase the 
likelihood of a merger of the two cores; in fact, by providing an upper bound on the infall 
rate, we see that both cores will form protostars long before merging. We have calculated 
the average density of the background medium using spheres both two and five times the 
radius of the core separation and found that the two cores should not approach closer 
than 100 AU for at least 3600 years. 

We plotted radially-binned, mass-weighted average quantities at several epochs of the 
halo collapse. (Fig. 3) The central bin was located at the most dense zone in the simu- 
lation at all times. After the cloud fragmented, this most dense zone was always within 
Core A. The location of Core B is clearly visible as a spike in the density, temperature, 
and H 2 fraction, and a slight enhancement in the enclosed mass, at ~ 800AU. The tem- 
perature plot shows a marked increase at high densites compared to that of both previous 
studies (1,11). This is likely a result not only of different molecular hydrogen formation 
and cooling rates, and also the inclusion of the heating from the formation of molecular 
hydrogen and opacity to molecular hydrogen lines. These factors delay the conversion to 
a fully-molecular state until the gas collapses to higher densities, reducing the efficiency 
of the ro- vibrational cooling, leading to higher temperatures. Furthermore, the average 
molecular hydrogen fraction of the cloud does not reach unity, as it has in previous calcu- 
lations. At higher densities the cooling from ro- vibrational lines is inefficient as a result of 
increased optical depth; this leads to an increase in the temperature, and thus a drop in 
the molecular hydrogen fraction. There is no substantial dissociation at the final output 
of this simulation, but the inner region of the collapsing halo has a molecular hydrogen 
mass fraction of ~ 0.5. As the cloud collapses further, H 2 will dissociate at densities 
just slightly higher than those presented here. Following the prescription of (9), the peak 
density of the simulation provides an effective transmission coefficient of 2 x 10~ 3 for ro- 
vibrational cooling and a transmission coefficient of 0.94 for cooling from collision induced 
emission. 

To estimate instabilities brought on by chemical and radiative mechanisms, we follow 
the methods of (9); when the characteristic timescales for the change in energy of the gas is 
less than the dynamical time, small perturbations can no longer be effectively suppressed 
and the gas is susceptible to fragmentation. The dynamical time is defined by 



By defining E as the total thermal energy of the gas, e ra d as the radiative cooling from 
ro- vibrational lines and collision-induced emission, and e c h em as the heating or cooling due 
to chemical changes in the molecular hydrogen content of the gas (4.48 eV per molecule), 
we can write the molecular hydrogen timescale as 




E 




4 



Comparing these two quantities (Fig. 4) we found that the ratio was approximately unity 
at most densities, but that at the densities at which the gas fragmented, this ratio changed 
sign abruptly. This is a steep transition, from being moderately prone to fragmentation, 
to being dominated by the heating from the formation of H 2 . However, as shown by (22), 
in isolation the thermal instability brought on by the onset of three-body heating and 
molecular line cooling would not be sufficient to cause the fragmentation of a collapsing 
gas cloud. 

We define the rotational energy as 

E rot = -U)]W 

where u is the rotational velocity and I is the inertial tensor of the gas. For a sphere of 
radius r and mass M the gravitational binding energy W is 

W S ?^ 
5 r 

We compared these two quantities (Fig. 4, bottom) and indicated with thin lines the criti- 
cal ratios of ~ 0.44 (23) and 0.27 (24) above which dynamical instabilites can become dom- 
inant in compressible toroids and compressible spheroids with a centrally-concentrated 
mass. This ratio was at a maximum at an enclosed mass of approximately 1OOM . The 
transition between instabilities identified above is within this radius. This halo fragments 
into multiple components because of the nesting of this marginal chemical instability 
within the radius of that of the rotational instability. However, more important is if an 
increase in this ratio is found exactly at the radii where the cloud has been identified as 
fragmenting. While analytic estimates for the gravitational potential energy of a cloud 
vary by factors of order unity depending on the eccentricity and ellipticity of the cloud, 
the relative ratio traces the mass scales where fragmentation may occur. 

The total number of first-generation Population III stars is likely to be significantly 
smaller than that of the second generation. Assuming no recombination, massive metal- 
free stars can ionize of order 10 7 M of gas, which is between 100 — 1000 times the gas mass 
of a Population III star-forming halo (25), and thus will be outnumbered by a large factor 
by metal-free stars forming in preprocessed regions. Fragmentation in unprocessed, first- 
generation Population III stars thus indicates that fragmentation in preprocessed, second- 
generation Population III stars is also possible. This would ultimately be more important 
to structure formation and the star formation history of the universe and the chemical 
signature of metal free stellar evolution in the fossil record of old stars in the Galaxy. (26) 
By changing the size of the gas reservoirs from which the first-generation Population III 
protostars stars accrete, as well as the ultraviolet flux from the subsequent metal-free 
stars themselves, the long-term star formation environment for the next generation of 
Population III stars could be dramatically changed. A higher flux of ionizing photons 
could lead to more efficient cooling by molecular hydrogen in neighboring halos (27, 28) 
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but a lower meta-galactic flux would likely lead to a higher global star formation rate and 
higher accretion rates. (29) 

We have shown that in a particular realization fragmentation occurs at relatively high 
densities. The separation and timescales are such that the two identified fragments are 
likely to form a binary stellar system. Improvements to chemical solvers at high densities 
and better laboratory values for the molecular hydrogen rate coefficients may change 
the details of the collapse, but these results demonstrate that fragmentation is possible 
in metal-free halos collapsing from realistic initial conditions. The frequency of such 
fragmentation, and thus that of metal-free binary systems, cannot be gauged using the 
small number of simulations conducted thus far. 

The problem of "finding" fragmentation in cosmological halos may be one of poor 
sampling; if fragmentation is rare, the small number of published calculations likely will 
not sample those halos in which it could occur. In particular, if fragmentation is more 
likely to occur in halos that undergo rapid merger history, the current means of ab initio 
simulation of Population III star formation may be ill-equipped to study its likelihood 
and relevance. We have found fragmentation in one out of five realizations, suggesting 
that binaries are a formation channel that must be considered in population synthesis 
studies. From a set of only five realizations an initial mass function cannot be estimated. 
However, we conclude by noting that the current means of conducting high-dynamic range 
simulations are biased against finding fragmentation in collapsing pre-stellar clouds. This 
simulation fragmented at a density of ~ 2.0 x lCT 13 g cur 3 , which corresponds to a 
dynamical time of 210 years; if the time delay between clump formation is greater than 
approximately the dynamical time, the current methods of simulating primordial star 
formation will not observe fragmentation. Consequently it is likely that a substantial 
fraction of Population III stars are forming binaries or event multiple systems. 
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Fig. 1. 

Mass- weighted average density (left column), H2 fraction (middle column), and temper- 
ature (right column) projected through a cube centered on the center of mass of the 
two-core system with a side length of 3500 AU. The bottom row is the final output of 
the simulation, the middle row is 555 years previous, and the top row is an additional 591 
years previous to the middle row (1146 years before the end of the simulation). Gravita- 
tionally bound cores with minimum density 3.0 x 10~ 12 g cm -3 have been outlined with 
thick lines in the bottom row; Core A is on the left and Core B is on the right. Field of 
view is 3500AU. 



7 




°0 20 40 60 80 100 120 140 160 180 

t [years] 



Fig. 2. 

Cores, identified by placing a minimum density cut at 3.0 x lCT 13 g cm -3 , plotted over 
time. The y-position denotes current maximum density, x-position denotes time since 
the first core was identified. Overlaid is the predicted free-fall evolution of the peak den- 
sity in each core starting from the first time when the m aximum density is greater than 
lCT 12 g cm -3 , governed by ^ = where we use tff = \J 32 Qp . We have extended the 
free-fall trajectories past the end of the simulation, which is indicated by the dotted line. 
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Fig. 3. 

Mass-weighted, spherically-averaged quantities as a function of distance from the most 
dense point (which is located within Core A after the cloud fragments) at different times 
in the simulation: mass-enclosed as a function of radius (Panel A), density as a function 
of radius (Panel B), molecular hydrogen mass fraction (Panel C), and temperature (Panel 
D). 
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Fig. 4. 

Enclosed quantities as a function of mass-enclosed, with respect to the most dense point 
(which is located within Core A after the cloud fragments) and calculated in the rest 
frame of that point at different times in the simulation: the mass-weighted average ratio 
of the dynamical time of the gas divided by the cooling time of the molecular hydrogen, 
taking into account the heating from three-body formation processes (top panel), and 
rotational energy divided by gravitational binding energy (bottom panel). In the bottom 
panel, lines have been drawn to indicate a ratio of 0.27 (thin, dashed) and 0.44 (thin, 
solid). 

Supporting Online Material 
Simulation Methodology 

The simulation presented in this work corresponds to simulation L0_30C in (11), which 
was initialized from Gaussian random initial conditions at redshift z = 99 inside a cubic 
volume with a side length of 300 kpc h _1 (comoving). In addition to standard refinement 
criteria based on overdensity and cooling time, we require a minimum of 64 cells per 
Jeans length, ensuring adequate resolution at all stages of collapse, ending with 31 levels 
of refinement. 

The simulations presented here implement a primordial chemistry solver optimized 
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to be efficient in the stiff regime when three-body reactions forming molecular hydrogen 
become dominant, utilizing the method from (30). The equation of state of molecular hy- 
drogen is also treated accurately, and all kinetic rate equations are solved without applying 
prescriptions for equilibria. In contrast to previous adaptive mesh refinement simulations, 
these calculations include energy contributions due to the formation and destruction of 
molecular hydrogen as a result of three-body reactions (4.48 eV/molecule), the dominant 
H 2 formation channel above densities of 10~ 16 g cm" 3 . These energy contributions, as well 
as energy losses due to radiative processes, are split into formation and destruction con- 
tributions and the energy is updated in a second-order backward difference method. This 
allows for fast numerical convergence during integration. The changes in the gas chem- 
ical state are updated on the same timescale as the thermal changes, ensuring accurate 
chemothermal composition. 

For the coefficients governing the formation and dissociation of molecular hydrogen, 

(k 22 ) H + 2H -> H + H 2 
(£43) H + H 2 -> H + 2H 

we use the rate calculated in (10) where a discussion of the uncertainty governing 
this particular process is detailed. This uncertainty is further discussed in (31). We 
have also used a prescription from (9) for the optical depth of ro- vibrational cooling from 
molecular hydrogen, as well as the optical depth to collision-induced absorption. We have 
additionally included the formation of molecular hydrogen via three-body reactions with 
H 2 as the third body, using the rates given in (32). The remaining set of chemical rate 
coefficients is identical to that of (11). 

We have simulated five different cosmological halos, with four of the calculations cor- 
responding to to runs from (11). The simulation presented in this work corresponds to 
simulation LCL30C in (11), which was initialized from Gaussian random initial conditions 
at redshift z = 99 inside a cubic volume with a side length of 300 kpc h -1 (comoving). The 
other four realizations did not fragment before we stopped each simulation. We initialize 
the simulation on a 128 3 root grid with three nested levels of refinement, for an effective 
resolution of 1024 3 (0.42 kpc proper at z — 99). Our finest dark matter particle mass is 
2.6M . The first protostar forms at redshift z = 19.08, inside a dark matter halo of mass 
5.8 x 10 5 M with a gas spin parameter of 0.025 and dark matter spin parameter of 0.042, 
where the spin parameter is defined as: 

_ L|£|V 2 
~ CM 3 / 2 

where L is the mass-weighted average specific angular momentum, E is the total kinetic 
energy and M is the enclosed mass. This collapse time is somewhat delayed from that 
in (11), which we attribute to different choices of molecular hydrogen formation rate and 
cooling rates; however, the two simulations are in general agreement. The simulation was 
stopped at a maximum (proper) baryon density of 1.61 x 10~ 8 g cm -3 («h ~ 10 16 cm~ 3 ), 
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at a redshift of 19.08, with 31 levels of refinement. The peak spatial resolution is 1.72 R & 
(proper) and at the end of the simulation there were 8.75 x 10 7 ~ 440 3 unique computa- 
tional elements with a final mass resolution of ~ 10 -7 M . 

Quantities 

We tabulate several relevant quantities of the two cores in Table [TJ These values have all 
been calculated during the final output from the simulation, when the maximum baryon 
density was 1.61 x 10~ 8 g cm -3 . 

Resolution Studies 

We have conducted resolution studies of this particular realization. With the refinement 
criteria such that we had 16 cells per Jeans length (J16), we allowed the simulation 
to run until the maximum density was 10~ 5 g cm -3 , where we observed that two cores 
had formed with a separation of ~ 200 AU and masses of two and seven solar masses, 
with topologically distinct, gravitationally-bound density contours at a density of 9.5 x 
10~ 13 g cm -3 . The separation was found to be lower than the Keplerian radius. To 
ensure that this fragmentation was not the result of resolution-limited fragmentation, 
we restarted the simulation from the time when the maximum density was 10~ 21 g cm -3 , 
before the onset of rapid molecular hydrogen formation. With a total of 32 cells per Jeans 
length (J32), we saw fragmentation at length scales of ~ 100 AU. However, because of 
the close pair distance, it is unclear whether or not the cores would subsequently merge. 
The work presented here had the greatest resolution, requiring 64 cells per Jeans length 
(J64). 

By changing the required resolution, we have changed the structure of the hydrody- 
namic flow into the collapsing region. The delay time in the simulation corresponded 
roughly to the strength of fragmentation; J32 collapsed the earliest, with J16 collapsing 
6000 years later and J64 collapsing 175,000 years following that. This corresponds to 
roughly 8% of the freefall time at the density at which we restarted the simulation. The 
actual collapse takes about 5 million years because of the balance of molecular cooling 
and heat input from turbulence on the large scales (~ pc) as hitherto seen and explained 
by (20). As we conduct the resolution study varying the spatial resolution by a factor 
4 one also slightly modifies the "initial conditions" since we start the calculation not 
shortly after recombination, where the original simulation was initialized, but rather at 
an intermediate time. The loitering phase relevant for primordial star formation otherwise 
will always be prone to the intermittency issues known so well in interstellar turbulence. 
Nevertheless, looking at the difference between the J16 and J32 runs the slight change in 
"initial conditions" only lead to a 6000 year difference. The timing clearly is the most 
significant differences in these runs. 
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Core Finding Algorithm 

We identified cores by using an optimized version of the core finder described in (33). 
Starting with a sphere of radius 2500AU in radius, centered at the most dense point 
in the simulation, we recursively identify topologically connected sets of cells of a given 
density minimum, incrementing the density minimum with every subset of the data that is 
extracted. All sets of cells that are not gravitationally bound are eliminated, as are objects 
with only a single subset identified. This allows us to identify the largest topologically 
connected sets of cells that are gravitationally bound. 

In this algorithm we compute the gravitational energy by direct summation, rather 
than using the potential Enzo computes. Consequently the computational cost scales 
with N 2 where N is the number of cells in the simulation. As such, performing this 
computation on commodity CPUs became cost prohibitive; to that end, we converted 
the code to run on the NVidia CUDA framework^ This decreased computational time 
substantially, allowing us to perform finely-spaced sweeps of density-space (separated by 
0.15 dex), ensuring we have accurately identified fragments. 
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Core p max Mass R tff 

~A 1.61 x 10~ 8 g cm- 3 1O.OM 24.4 AU 0.52 years 
B 5.03 x 10- 12 g cm- 3 6.3M 98 AU 29.7 years 

Table 1: Table of quantities for identified gas cores from final simulation output 
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